@Gabriele Ghisleni, Data Science.
Install the required dependecies on google colab and mount the drive.
%%capture
!pip install geopandas==0.10.2
!pip install mapclassify==2.4.3
!pip install pygeos==0.10.2
!pip install rasterio==1.2.3
!pip install owslib==0.25.0
!pip install -U folium
import warnings, sys, os, folium, owslib, rasterio, geopy, shapely, json, requests, pprint
from folium import FeatureGroup, GeoJson, Marker
import matplotlib.pyplot as plt
from google.colab import drive
from tqdm.auto import tqdm
import geopandas as gpd
import pandas as pd
import numpy as np
drive.mount('/content/drive')
warnings.simplefilter("ignore")
p_print = pprint.PrettyPrinter(indent=4).pprint
# change this to your folder directory
folder_directory = "/content/drive/MyDrive/_____SHARED/Geospatial"
sys.path.append(folder_directory)
%%capture
#download the necessary utils files
if 'geo_utils' not in os.listdir():
os.mkdir('geo_utils')
if "utils.py" not in os.listdir('geo_utils'):
url_utils = "https://raw.githubusercontent.com/GabrieleGhisleni/Geo-Dolomities/master/utils/utils.py"
url_style = "https://raw.githubusercontent.com/GabrieleGhisleni/Geo-Dolomities/master/utils/styles.py"
!wget --no-cache --backups=1 -P'utils' {url_utils}
!wget --no-cache --backups=1 -P'utils' {url_style}
!touch utils/__init__.py
assert 'utils.py' in os.listdir('geo_utils'), 'utils.py not present'
assert 'styles.py' in os.listdir('geo_utils'), 'styles.py not present'
#in case ModuleError restart the kernel
import geo_utils.utils as utils
import geo_utils.styles as styles
tiles = 'Stamen Terrain'
drive.mount('/content/drive')
folder_directory = "/content/drive/MyDrive/_____SHARED/Geospatial"
sys.path.append(folder_directory)
os.chdir(folder_directory)
import utils, styles
tiles = 'Stamen Terrain'
Mounted at /content/drive
Let's explore the area of the dolomities, the data are from the offical website of Dolomiti Unesco, that are available in .kml format.
We load the data and perform some preprocessing before so to have a better insight.
# read kml in GeoPandas https://docs.astraea.earth/hc/en-us/articles/360043923831-Read-a-KML-File-into-a-GeoPandas-DataFrame
gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'
dolomiti_df = gpd.read_file('data/geo_dolomiti/dolomiti_area.kml', driver='KML')
print(f"Reference system information:\n{ '-'*50}\n{dolomiti_df.crs.__repr__()}{'-'*50}")
Reference system information: -------------------------------------------------- <Geographic 2D CRS: EPSG:4326> Name: WGS 84 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: World. - bounds: (-180.0, -90.0, 180.0, 90.0) Datum: World Geodetic System 1984 ensemble - Ellipsoid: WGS 84 - Prime Meridian: Greenwich --------------------------------------------------
# some preprocessing
dolomiti_df["name"] = dolomiti_df["Name"].apply(lambda x: x.replace(',', ' &').replace('\xa0', ' '))
dolomiti_df.loc[2, 'name'] = "Sistema 3 - Pale di San Martino"
dolomiti_df['url_info'] = dolomiti_df.Description.apply(lambda x: x.split('>')[1])
# representative points and area expressed in m^2 (project to 32632 which is in meters).
dolomiti_df['representative_point'] = dolomiti_df.to_crs(epsg=32632).geometry.representative_point().to_crs(epsg=4326)
dolomiti_df['representative_point'] = dolomiti_df['representative_point'].apply(str)
dolomiti_df['area'] = dolomiti_df.to_crs(epsg=32632).geometry.area.apply(lambda x: int(x)) #metersquare since we are in 32632
dolomiti_df = dolomiti_df.loc[:, ['name', 'url_info', 'area', 'representative_point', 'geometry']]
dolomiti_df.head(3)
| name | url_info | area | representative_point | geometry | |
|---|---|---|---|---|---|
| 0 | Sistema 1 - Pelmo & Croda da Lago | http://www.dolomitiunesco.info/?gruppo-dolomit... | 43153031 | POINT (12.08268734287994 46.45014333277591) | POLYGON Z ((12.12142 46.40488 0.00000, 12.1214... |
| 1 | Sistema 2 - Marmolada | http://www.dolomitiunesco.info/?gruppo-dolomit... | 22057860 | POINT (11.859039731003216 46.428646792973986) | POLYGON Z ((11.82319 46.45296 0.00000, 11.8221... |
| 2 | Sistema 3 - Pale di San Martino | http://www.dolomitiunesco.info/?gruppo-dolomit... | 316567741 | POINT (12.165681166978791 46.23429736260587) | POLYGON Z ((11.98683 46.31720 0.00000, 11.9880... |
# Some information regarding the groups and their extensions
for idx, row in dolomiti_df.sort_values('area', ascending=False).reset_index().iterrows():
print(f"""\
{idx+1}° largest group: '{row["name"]}' Area: {round(row.area * 1e-6)} km2.
Website for more info: {row.url_info}.
""")
1° largest group: 'Sistema 5 - Dolomiti settentrionali' Area: 536 km2. Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=dolomiti-settentrionali. 2° largest group: 'Sistema 3 - Pale di San Martino' Area: 317 km2. Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=pale-di-san-martino. 3° largest group: 'Sistema 4 - Dolomiti friulane e d'Oltre Piave' Area: 215 km2. Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=dolomiti-friulane. 4° largest group: 'Sistema 9 - Dolomiti di Brenta' Area: 111 km2. Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=dolomiti-di-brenta. 5° largest group: 'Sistema 7 - Sciliar-Catinaccio & Latemar' Area: 93 km2. Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=sciliar-catinaccio-latemar. 6° largest group: 'Sistema 6 - Puez-Odle' Area: 79 km2. Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=puez-odle. 7° largest group: 'Sistema 1 - Pelmo & Croda da Lago' Area: 43 km2. Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=pelmo-croda-da-lago. 8° largest group: 'Sistema 2 - Marmolada' Area: 22 km2. Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=marmolada. 9° largest group: 'Sistema 8 - Bletterbach' Area: 3 km2. Website for more info: http://www.dolomitiunesco.info/?gruppo-dolomitico=bletterbach.
# first representation and plot of the dolomities
fig, axes = plt.subplots(1,1, figsize=(15,15))
dolomiti_df.to_crs(epsg=4326).plot(
ax=axes,
column='name',
cmap = "tab10",
legend=True,
categorical=True,
legend_kwds= {
'bbox_to_anchor':(.3, 1.05),
'fontsize':15,
'frameon':False
}
)
plt.title('Dolomities Area', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()
Extract the regions of the italy that are closeby the dolomities.
italy_regions = gpd.read_file("data/istat_administrative.gpkg", layer="regions")
# Set correct crs
italy_regions = italy_regions.to_crs("EPSG:4326")
# Selecting only interested regions
assert italy_regions.crs == dolomiti_df.crs
selected_regions = []
for _, row in dolomiti_df.iterrows():
regions = italy_regions.geometry.overlaps(row.geometry)
selected_regions.extend(
italy_regions.loc[regions].DEN_REG.values.tolist()
)
selected_regions = list(set(selected_regions))
dolomiti_regions = italy_regions.loc[italy_regions.DEN_REG.isin(selected_regions)]
Using Folium library for creating nice interactive maps so to explore the dolomities with respect to the italian regions.
# https://gis.stackexchange.com/questions/392531/modify-geojson-tooltip-format-when-using-folium-to-produce-a-map
# https://stackoverflow.com/questions/58305337/plotting-polygons-in-python-using-geopandas-and-folium
# https://stackoverflow.com/questions/63960271/folium-geojson-custom-color-map
# https://leafletjs.com/reference.html#path-option
dolomiti_map = utils.get_new_map(title='Dolomities in Italy', tiles=tiles)
feature_layers = dict(
dolomiti_area = FeatureGroup(name='Area Dolomiti'),
regions_layer = FeatureGroup(name='Italian Regions')
)
dolomiti_region = GeoJson(
data = dolomiti_regions,
style_function = styles.style_region,
tooltip = folium.GeoJsonTooltip(fields=['DEN_REG'])
).add_to(feature_layers["regions_layer"])
dolomiti_geo = GeoJson(
data = dolomiti_df,
style_function = styles.style_dolomiti,
tooltip = folium.GeoJsonTooltip(
fields=[
'name',
'url_info',
'area',
'representative_point'
])
)
dolomiti_geo.add_to(feature_layers["dolomiti_area"])
utils.add_map_infos(feature_layers, dolomiti_map, reverse_order=True)
# explore with static maps
fig, axes = plt.subplots(1,1, figsize=(15,15))
dolomiti_regions.to_crs(epsg=4326).plot(
ax=axes,
cmap = "Pastel1",
column='DEN_REG',
categorical=True,
legend=True,
)
dolomiti_df.to_crs(epsg=4326).plot(
ax=axes,
column='name',
cmap = "tab10",
legend=True,
categorical=True,
legend_kwds= {
'bbox_to_anchor':(.3, 1.05),
'fontsize':10,
'frameon':False
}
)
plt.title('Dolomities Area', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()
italy_municipalities = gpd.read_file("data/istat_administrative.gpkg", layer="municipalities")
# filter by the interested regions
italy_municipalities = italy_municipalities.loc[italy_municipalities.COD_REG.isin(dolomiti_regions.COD_REG.unique())]
# Set correct crs
italy_municipalities = italy_municipalities.to_crs("EPSG:4326")
# find the municipalities that are nearby the dolomities area
def geometry_in_areas(geom, regions_df, attribute_name='name'):
res = []
for idx, region in regions_df.iterrows():
geometry = region.geometry
if geom.overlaps(geometry):
res.append(region[attribute_name])
return ','.join(res)
buffered_dolomiti = dolomiti_df.copy()
buffered_dolomiti.geometry = dolomiti_df.to_crs(epsg=32632).geometry.buffer(500).to_crs(epsg=4326).values
selected_municipalities = italy_municipalities.geometry.apply(lambda x: geometry_in_areas(x, buffered_dolomiti)) # dolomiti_df dolomiti_regions
italy_municipalities['dolomiti'] = selected_municipalities.values
dolomiti_municipalities = italy_municipalities.loc[selected_municipalities != '']
dolomiti_municipalities = dolomiti_municipalities.loc[:, ['COD_REG', 'COMUNE', 'geometry', 'dolomiti']]
dolomiti_municipalities.head(2)
| COD_REG | COMUNE | geometry | dolomiti | |
|---|---|---|---|---|
| 2760 | 4 | Badia | MULTIPOLYGON (((11.89993 46.63847, 11.89978 46... | Sistema 5 - Dolomiti settentrionali,Sistema 6 ... |
| 2763 | 4 | Braies | MULTIPOLYGON (((12.13557 46.73782, 12.13594 46... | Sistema 5 - Dolomiti settentrionali |
# add regions of the municipality
dolomiti_municipalities = dolomiti_municipalities.merge(
dolomiti_regions.loc[:, ['COD_REG', 'DEN_REG']],
on='COD_REG'
)
# municipality for each region that touches the dolomities
for _, row in dolomiti_municipalities.groupby('DEN_REG').size().sort_values(0, ascending=False).reset_index().iterrows():
print(f"'{row.DEN_REG}' has {row[0]} municipalities nearby the dolomities.")
'Trentino-Alto Adige' has 44 municipalities nearby the dolomities. 'Veneto' has 40 municipalities nearby the dolomities. 'Friuli Venezia Giulia' has 6 municipalities nearby the dolomities.
# create the map and add geojson infos
dolomiti_map_municipalities = utils.get_new_map(title='Dolomities over municipalities', tiles='OpenStreetMap')
feature_layers = dict(
dolomiti_area = FeatureGroup(name='Area Dolomiti'),
municipalities_layer = FeatureGroup(name='Italian Municipalities')
)
dolomiti_geo.add_to(feature_layers["dolomiti_area"])
GeoJson(
data = dolomiti_municipalities,
style_function = styles.style_municipalities,
tooltip = folium.GeoJsonTooltip(
fields=[
'COMUNE',
'DEN_REG',
'dolomiti'
])
).add_to(feature_layers["municipalities_layer"])
utils.add_map_infos(feature_layers, dolomiti_map_municipalities, reverse_order=True)
Manually extract the CAI from the CAI official website alpine hut and organize them into a useful way.
# https://www.cai.it/andare-in-montagna/rifugi-e-bivacchi/rifugi-del-club-alpino-italiano/
rifugi_trentino = pd.read_excel(
io='data/alpine_huts/rifugi_trentino.xlsx',
sheet_name=0,
header=1
)
rifugi_veneto = pd.read_excel(
io='data/alpine_huts/rifugi_veneto.xlsx',
sheet_name=0,
header=1
)
rifugi = pd.concat([rifugi_trentino, rifugi_veneto])
rifugi.head(2)
| Nome | Tipo | Sezione | Quota (m slm) | Posti totali | Categoria | |
|---|---|---|---|---|---|---|
| 0 | Rifugio Vioz | Rifugio custodito | S.A.T. | 3535 | 66 | E |
| 1 | Rifugio Boè | Rifugio custodito | S.A.T. | 2873 | 69 | D |
user_agent = "Mozilla/5.0 (Linux; Android 10) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/86.0.4240.75 Mobile Safari/537.36"
geolocator = geopy.geocoders.Nominatim(user_agent=user_agent, timeout=3)
def geocode_sat(rifugi):
detail = []
not_founded = []
for idx, row in tqdm(rifugi.iterrows(), total=len(rifugi)):
name = row['Nome']
# manually fix the names that differ
if name in utils.names_adj:
name = utils.names_adj[name]
# geolocate the huts
geocodes = geolocator.geocode(
query = name,
addressdetails = True,
extratags = True
)
if geocodes:
detail.append(geocodes.raw)
else:
not_founded.append(row['Nome'])
rifugi = utils.expand_raw_geocode(pd.DataFrame(detail))
rifugi = gpd.GeoDataFrame(rifugi, crs="EPSG:4326")
rifugi = rifugi.loc[~rifugi.county.isin(['Roma Capitale', "Reggio nell'Emilia", 'Bergamo'])]
rifugi.drop(columns=['boundingbox']).to_file('data/rifugi.shp')
return rifugi
cai_rifugi = (
gpd.read_file("data/alpine_huts/rifugi.shp")
if os.path.exists('data/alpine_huts/rifugi.shp')
else geocode_sat(rifugi)
)
cai_rifugi.head(2)
| municipali | county | state | road_name | extended_n | opening_ho | contact_mo | winter_roo | capacity | website | display_na | type_descr | geometry | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Comunità della Valle di Sole | Provincia di Trento | Trentino-Alto Adige/Südtirol | Vioz | Rifugio Mantova al Vioz | Jun 20-Sep 20 | +39 339 2798826 | yes | 30 | http://www.rifugiovioz.it/ | info@rifugiovioz.it | Rifugio Mantova al Vioz, Vioz, Peio, Comunità ... | alpine_hut | POINT (10.63578 46.39920) |
| 1 | Pustertal - Val Pusteria | Bolzano - Bozen | Trentino-Alto Adige/Südtirol | Cresta Strenta - Lichtenfelser Weg | Rifugio Boé | Jun 20-Sep 20 | +39 349 5730110 | yes | 69 | http://www.rifugioboe.it/ | rifugioboe@gmail.com | Rifugio Boé, Cresta Strenta - Lichtenfelser We... | alpine_hut | POINT (11.82328 46.51461) |
sat_map = utils.get_new_map(title='Hut mantained by SAT', tiles=tiles)
unique_provice_dolomities = cai_rifugi.county.unique().tolist()
colors = ['gray', 'lightgreen', 'green', 'orange', 'darkred','black']
feature_layers = {
group:FeatureGroup(name=f"<span style='color:{utils.c_colors[idx]}'>Hut in {group}</span>")
for idx, group in enumerate(unique_provice_dolomities)
}
feature_layers["dolomiti_area"] = FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers["dolomiti_area"])
for idx, row in cai_rifugi.iterrows():
marker = Marker(
location = [i[0] for i in reversed(row.geometry.xy)],
popup = utils.get_info_from_row(row),
tooltip = utils.get_info_from_row(row),
icon = folium.Icon(color= utils.c_colors[unique_provice_dolomities.index(row.county)])
)
marker.add_to(feature_layers[row.county])
utils.add_map_infos(feature_layers, sat_map)
# explore with static maps
fig, axes = plt.subplots(1,1, figsize=(15,15))
dolomiti_regions.to_crs(epsg=4326).plot(ax=axes, linewidth=3, facecolor='none', edgecolor='gray')
dolomiti_df.to_crs(epsg=4326).plot(
ax=axes,
linewidth=2,
edgecolor='black',
column='name',
cmap = "tab10",
legend=True,
categorical=True,
legend_kwds= {
'bbox_to_anchor':(.3, 1.05),
'fontsize':10,
'frameon':False
}
)
cai_rifugi.to_crs(epsg=4326).plot(ax=axes, color = "orange")
plt.title('Huts over the Dolomities area', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()
from owslib.csw import CatalogueServiceWeb
# fetch the catalog
csw = CatalogueServiceWeb("http://www.pcn.minambiente.it/geoportal/csw")
query = owslib.fes.PropertyIsLike('csw:AnyText', 'Laghi')
csw.getrecords2(constraints=[query], maxrecords=1_000)
#explore the available data
for i, rec in enumerate(csw.records):
print(f"{rec} - {csw.records[rec].title}")
if i > 15: break
r_campan:16050024:20090806:141714 - Specchi acqua campania r_piemon:4ba5287d-0c81-4729-84b9-906ad5ccb173 - Idrografia 1:100.000 - Laghi - Storico r_vda:03151-META:20170302:081500 - Arpa - Catasto laghi della Valle d'Aosta r_piemon:60048c12-904a-4aa5-a062-1a9dd94ae971 - Ppr - Laghi (tavv. P4-P5) p_TN:0403efca-cadb-4463-a656-4d0fba207953 - Aree di protezione dei laghi - Variante PUP 2000 r_piemon:4fb63c5b-813b-4c6d-9fed-6ce92df6ca64 - Ppr - Laghi (tav. P2) r_basili:54eecefc:162a68d9be7:-2723 - Beni paesaggistici art. 142 let. b - Laghi ed invasi artificiali (Buffer) ispra_rm:20101021:113000 - Reticolo Idrografico Nazionale - Dataset r_lombar:36BB78E0-8946-4E7D-ADAC-7D149D18B570 - PTUA 2016 - Tav. 1 - Corpi idrici superficiali e bacini drenanti (fiumi e laghi) R_SARDEG:XBDOJ - PPR06 - Laghi naturali, invasi artificiali, stagni e lagune r_campan:16050003:20090722:144054 - Aree di tutela coste e laghi r_piemon:b65fc959-dfa5-4edd-991d-b6b9f638ab5b - Ppr - Lettera b) fasce lacuali (tav. P2) p_TN:9a5dc19c-543f-44d3-a36c-b4ea73d0a48a - Z101_P_PUP abnfa_tn:00002:20120504:130718 - Laghi Alpi Orientali p_TN:822c7446-8e84-45b6-9d50-98d5dca19416 - Z310_P_PUP r_sicili:aa3781a2-3c7d-4a32-ac9d-e8aae4195f1f - Carta dell' uso dei suoli 1:100.000 r_piemon:8d267850-e543-4294-bb23-483aa5335183 - Corpi idrici WFD 2000/60/CE Triennio 2009-2011 - Geo servizio WFS
Explore the lakes in the dolomities.
# explore the resource found
unique_id = "m_amte:299FN3:5237cf1e-4d49-486c-a688-152bc4473508"
dicts = csw.records[unique_id].__dict__
print('Abstract:\n'); print(dicts['abstract']); print('-'*100)
print('References:\n'); p_print(dicts['references']); print('-'*100)
print('Information:\n'); p_print({f"{k}: {v}" for k,v in dicts.items() if v and k != 'xml' and k != 'abstract' and type(v) != list}); print('-'*100)
Abstract: Elementi poligonali appartenenti al DBPrior10K rappresentanti gli specchi dacqua interni come laghi, lagune, stagni, bacini artificiali, ecc. ---------------------------------------------------------------------------------------------------- References: ---------------------------------------------------------------------------------------------------- Information: ----------------------------------------------------------------------------------------------------
from owslib.wfs import WebFeatureService
url = "http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/Vettoriali/Specchi_Acqua.map"
wfs = WebFeatureService(url=url, version="1.1.0")
{f"{k}: {v}" for k,v in wfs.__dict__.items() if v and k != 'xml' and k != 'abstract' and type(v) != list}
{'_capabilities: <Element {http://www.opengis.net/wfs}WFS_Capabilities at 0x7f061f7c6460>',
'auth: <Authentication shared=False username=None password=None cert=None verify=True auth_delegate=None>',
"contents: {'ID.ACQUEFISICHE.SPECCHI.ACQUA': <owslib.feature.wfs110.ContentMetadata object at 0x7f061e43e9d0>}",
'filters: <owslib.fes.FilterCapabilities object at 0x7f061d908b90>',
'identification: Service: OGC WFS, title=',
'owscommon: <owslib.ows.OwsCommon object at 0x7f061d90e250>',
'provider: <owslib.ows.ServiceProvider object at 0x7f061d9846d0>',
'timeout: 30',
'url: http://wms.pcn.minambiente.it/ogc?map=%2Fms_ogc%2FWMS_v1.3%2FVettoriali%2FSpecchi_Acqua.map',
'version: 1.1.0'}
# load the layer as geodataframe
lakes_df = gpd.read_file(
wfs.getfeature(
typename="ID.ACQUEFISICHE.SPECCHI.ACQUA",
bbox = dolomiti_df.total_bounds.tolist(),
)
)
# invert the xy coordinates
def reverse_geom(geom):
def _reverse(x, y):
return y[::-1], x[::-1]
return shapely.ops.transform(_reverse, geom)
lakes_df.geometry = lakes_df.geometry.apply(reverse_geom)
buffered_dolomiti = dolomiti_df.copy()
buffered_dolomiti.geometry = dolomiti_df.to_crs(epsg=32632).geometry.buffer(5_000).to_crs(epsg=4326).values
# find the lakes in the dolomities area
idxs = []
for idx, row in lakes_df.to_crs(epsg=4326).iterrows():
if buffered_dolomiti.geometry.contains(row.geometry).sum() > 0:
idxs.append(idx)
lakes_df = lakes_df.iloc[idxs]
lakes_df = lakes_df.loc[lakes_df.nome != ""].reset_index(drop=True)
lakes_df["area"] = lakes_df.to_crs(epsg=32632).geometry.area
lakes_df["area"] = lakes_df.to_crs(epsg=32632).geometry.area
lakes_df["area"] = (lakes_df["area"] * 1e-6).apply(lambda x: round(x, 3))
#biggest lakes
for _, row in lakes_df.sort_values('area', ascending=False).head().iterrows():
print(f"Name: {row['nome']} area: {row['area']}")
Name: LAGO DI MOLVENO area: 3.218 Name: Lago del Mis area: 1.098 Name: LAGO ARTIFICIALE DI FEDAIA area: 0.545 Name: LAGO DI TOVEL area: 0.371 Name: LAGO DELLO SCHENER area: 0.356
Exploration of all the glaciers that are in the dolomities.
data from http://www.glaciologia.it/en/, thanks to:
Salvatore M.C., Zanoner T., Baroni C., Carton A., Banchieri F.A., Viani C., Giardino M., Perotti L. (2015) - The state of Italian glaciers: A snapshot of the 2006-2007 hydrological period. Geografia Fisica e Dinamica Quaternaria, 38 (2), pp. 175-198.
# load the data if present otherwhise download
if 'ghiacciai.json' not in os.listdir('data'):
raw_glacier = requests.get('https://repo2.igg.cnr.it/ghiacciaiCGI/data/datagralp.json').json()
else:
with open('data/ghiacciai.json') as f:
raw_glacier = json.load(f)
# set the correct crs
crs = raw_glacier['crs']['properties']['name'].split(':')[-1]
glacier_df = gpd.GeoDataFrame.from_features(raw_glacier, crs=crs)
# keeps only the glaciers into the dolomities
idxs = []
for idx, row in glacier_df.to_crs(epsg=4326).iterrows():
if dolomiti_df.geometry.contains(row.geometry).sum() > 0:
idxs.append(idx)
glacier_df = glacier_df.iloc[idxs].reset_index(drop=True)
glacier_df["area"] = glacier_df.to_crs(epsg=32632).geometry.area
glacier_df["area"] = (glacier_df["area"] * 1e-6).apply(lambda x: round(x, 3))
#biggest glacier
for _, row in glacier_df.sort_values('area', ascending=False).head().iterrows():
print(f"Name: {row['GLACIER']} area: {row['area']}, gruppo: {row['GROUP_NAME']}")
Name: Ghiacciaio Principale della Marmolada area: 1.48, gruppo: Marmolada Name: Ghiacciaio Superiore (Orientale) dell'Antelao area: 0.255, gruppo: Antelao Name: Ghiacciaio (Vedretta) del Lagol o de l'Agola area: 0.215, gruppo: Brenta Name: Ghiacciaio di Travignolo area: 0.194, gruppo: Pale di San Martino Name: Ghiacciaio Occidentale del Sorapis area: 0.179, gruppo: Sorapiss
areas_glacier = glacier_df.groupby('GROUP_NAME').sum().reset_index().loc[:, ['GROUP_NAME', 'area']].sort_values('area', ascending=False)
for _, row in areas_glacier.iterrows():
print(f'Glacier group: {row.GROUP_NAME} has {round(row.area, 3)} km^2 of glaciers ')
Glacier group: Marmolada has 2.025 km^2 of glaciers Glacier group: Brenta has 1.124 km^2 of glaciers Glacier group: Pale di San Martino has 0.602 km^2 of glaciers Glacier group: Sorapiss has 0.453 km^2 of glaciers Glacier group: Antelao has 0.404 km^2 of glaciers Glacier group: Marmarole has 0.384 km^2 of glaciers Glacier group: Cristallo has 0.335 km^2 of glaciers Glacier group: Popera has 0.263 km^2 of glaciers Glacier group: Pelmo has 0.176 km^2 of glaciers Glacier group: Civetta has 0.109 km^2 of glaciers Glacier group: Tofane has 0.062 km^2 of glaciers Glacier group: Fanis has 0.037 km^2 of glaciers Glacier group: Croda Rossa d'Ampezzo has 0.001 km^2 of glaciers Glacier group: Moiazza has 0.001 km^2 of glaciers
# explore the lakes and the glaciers in the dolomities
glacier_lake_map = utils.get_new_map(title='Glaciers and Lakes', tiles='Stamen Terrain')
feature_layers = dict(
dolomiti_area = FeatureGroup(name='Area Dolomiti'),
glacier = FeatureGroup(name='Glacier'),
lakes = FeatureGroup(name='Lakes'),
)
dolomiti_geo.add_to(feature_layers["dolomiti_area"])
GeoJson(
data = glacier_df,
style_function = lambda x : {'color': 'purple', 'fill': True},
tooltip = folium.GeoJsonTooltip(
fields=[
'GLACIER',
'GROUP_NAME',
'MAX_ELEV',
'MIN_ELEV',
'DATE',
'AREA_KM'
])
).add_to(feature_layers["glacier"])
GeoJson(
data = lakes_df,
style_function = lambda x : {'color': 'blue', 'fill': True},
tooltip = folium.GeoJsonTooltip(
fields=[
'nome',
'area'
])
).add_to(feature_layers["lakes"])
utils.add_map_infos(feature_layers, glacier_lake_map)
# explore the lakes and the glaciers in the dolomities
fig, axes = plt.subplots(1,1, figsize=(15,15))
dolomiti_df.to_crs(epsg=4326).plot(
ax=axes,
linewidth=2,
edgecolor='black',
column='name',
cmap = "tab10",
alpha=.15,
legend=True,
categorical=True,
legend_kwds= {
'bbox_to_anchor':(.3, 1.05),
'fontsize':10,
'frameon':False
}
)
glacier_df.to_crs(epsg=4326).plot(ax=axes, color = "purple")
lakes_df.to_crs(epsg=4326).plot(ax=axes, color = "blue")
plt.title('Glaciers and Lakes in the Dolomities area', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()
Quickly exploration of the ortophoto of trentino only.
La Provincia Autonoma di Trento ha la piena proprietà dell'ortofoto. Attribuzione 3.0 Italia (CC BY 3.0 IT).
# fetch the catalog
csw = CatalogueServiceWeb("http://geodati.gov.it/RNDT/csw")
trento_query = owslib.fes.PropertyIsLike('csw:AnyText', 'Trento')
csw.getrecords2(constraints=[trento_query], maxrecords=1_000)
#explore the available data
for rec in csw.records:
if 'orto' in csw.records[rec].title.lower():
print(f"{rec} - {csw.records[rec].title}")
agea:00129:20090724:090446 - Ortofotocarta Trento 2003 agea:00377:20090911:093144 - Ortofotocarta Trento 2008 agea:00128:20090724:085449 - Ortofotocarta Trento 1997 c_l378:ortofoto2009 - Ortofoto 2009 p_TN:fbbc1e07-0b8e-46c9-b961-a02d8bebb217 - Ortofoto PAT 1973 in scala di grigi p_TN:013ef530-ee77-49d2-8f95-035b27ab1f0a - Ortofoto Val di Sella 2007 RGB p_TN:44cd8a88-ba6e-4299-8075-df9996401006 - Ortofoto Padergnone 2007 RGB p_TN:e506ab2b-5eb6-452a-8ebb-102ffc4454b1 - Ortofoto Paneveggio 2007 RGB p_TN:054e1af3-2dd8-4496-8b83-261cccc2674c - Ortofoto Ravina 2007 RGB p_TN:f2e88f1b-05d9-4942-93ee-857a0a9e1f0b - Ortofoto PAT 2015 RGB
from owslib.wms import WebMapService
from rasterio.plot import show
orto_photo_path = "p_TN:f2e88f1b-05d9-4942-93ee-857a0a9e1f0b"
wms_trento = csw.records[orto_photo_path].__dict__['references'][3]['url']
wms = WebMapService(wms_trento)
trentino_map = wms.getmap(
layers=['ecw-rgb-2015'],
srs='EPSG:25832',
format='image/jpeg',
bbox=(609500.0, 5056790.0, 731610.0, 5158990.0),
size=(3500,3500)
)
fig, axes = plt.subplots(1,1, figsize=(15,15))
fig = show(
rasterio.MemoryFile(trentino_map).open(),
ax=axes,
title='Trentino Natural View'
)
plt.xticks([]); plt.yticks([])
plt.show()